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Without invoking the Markov approximation, we derive formulas for vibrational energy relaxation 
(VER) and dephasing for an anharmonic system oscillator using a time- dependent perturbation 
theory. The system-bath Hamiltonian contains more than the third order coupling terms since 
we take a normal mode picture as a zeroth order approximation. When we invoke the Markov 
approximation, our theory reduces to the Maradudin-Fein formula which is used to describe VER 
I properties of glass and proteins. When the system anharmonicity and the renormalization effect due 

to the environment vanishes, our formulas reduce to those derived by Mikami and Okazaki invoking 
the path-integral influence functional method [J. Chem. Phys. 121, 10052 (2004)]. We apply our 
f~| ' formulas to VER of the amide I mode of a small amino-acide like molecule, A''-methylacetamide, in 

, heavy water. 
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I. INTRODUCTION 



PR 

O , Vibrational energy relaxation (VER) and dephasing are fundamental properties of molecular dynamics, energy 
' transfer, and reactivity. Many experimental and theoretical studies have explored these fundamental processes in gas 
I phase, the liquid state, and in glasses and biomolecular systems Q. Though our methodology can be applied to any 
I ^ molecular system, we are primarily interested in addressing VER and dephasing in peptides or proteins. While recent 
] advanced experimental techniques using absorption spectra or time-resolved spectra can deduce the structure and 
T-H . dynamics of such a peptide or protein system, theoretical approaches are needed to clarify the mechanisms of VER 
' and dephasing underlying the experimental data. 

, The most standard approach to this problem is through the perturbation theory of quantum mechanics as initiated 
^0 ■ by Oxtoby 01 ■ Recently Hynes's group (3| and Skinner's group 4] thoroughly studied the VER and dephasing 
I properties of water (their target mode was the OH bond of HOD in heavy water) using this strategy. This approach 
is applicable to peptides or proteins as was first illustrated by Straub and coworkers . Derived from this strategy is 
\^ [ the use of the Maradudin-Fein formula (or its equivalent) , which was pursued by Leitner ^ and Straub and coworkers 
■ • This formula requires the normal modes of the system and the cubic anharmonic coefficients between the normal 
' modes. This methodology can provide a reasonable account of VER properties of peptides or proteins, but there are 
• I— I . several deficiencies: the most serious one is that it assumes the Markov properties of the system, so it cannot describe 
• the short time dynamics @ . Another problem is the determination of the "lifetime" width parameter P, Q, ■ We 
also want to describe the dephasing properties of the system^rucial to the interpretation of the experimental results; 
' it cannot be directly described by the MF formula (but see 

To meet these goals, we derive the formulas for VER and dephasing without assuming the Markov properties, i.e., 
without taking an infinite-time limit. As a result, we can avoid the annoying "width parameter" problem inherent 
5-H ' to the MF approach. In this sense, Mikami and Okazaki 10] took a similar path using the path-integral influence 
functional theory. We use a simple time-dependent perturbation theory of quantum mechanics, and derive the VER 
and dephasing formulas more easily. We find there is a difference between our formulas and theirs in terms of 
renormalization of the system Hamiltonian. Another difference is that our system oscillator is taken to be a cubic 
anharmonic oscillator, whereas their mode is a harmonic oscillator. This can affect the result when the formulas are 
applied to real systems with strong anharmonicity. 

This paper is organized as follows: In Sec. ^ we derive the VER and dephasing formulas for an anharmonic 
oscillator (mode) without assuming the Markov properties. In Sec. IIIII we apply our formulas to the amide I mode of 
A^-methylacetamide in heavy water, and discuss the numerical results and the limitations of our strategy. In Sec. lIVI 
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we summarize the paper. Several system parameters and coefficients in our formulas are defined in the Appendix. 

II. DERIVATION OF THE FORMULAS FOR VER AND DEPHASING 
A. System, Bath, and Coupling 

We take our Hamiltonian of a solvated peptide or protein to be 

n = Ho + v = Hs + HB + v = Hf +nf + nB + v, (i) 

V = -qsi:F-{T))+ql{g-{g)) = -qsST + qlSg, (2) 
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- 0-5^1 + ^, (6) 
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qs = qs- -rr ^qs-b. (7) 

Tis — Ti.^s '' +^/ is the renormalized system Hamiltonian representing a vibrational mode qs with cubic anharmonicity 
/, Hb the bath Hamiltonian representing solvent or environmental degrees of freedom with harmonic frequencies uja, 
and V the interaction Hamiltonian describing the coupling between the system and the bath. We have assumed that 
the interation can be Taylor expanded, and we have only included up to the second order in qs- Note that we need 
to renormalize the system to assure that (V) = where the bracket denotes the bath average throughout this paper. 
(For the definition of SJ- and SG, see Appendix 1X1) This is automatically satisfied in the case of bilinear coupling like 
the Caldeira-Leggett-Zwanzig model llj, but this is not usually the case. The system variable becomes qs instead 
of qs, and the system frequency does lus instead of 0^5. This is similar to previous treatments of the system-bath 
interaction in the literature 
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B. Perturbation theory for VER and dephasing 



Starting from the interaction picture of the von Neumann equation, we can expand the density operator for the full 
system as 

m = pio) + \ f dt'[v{t'),~p{t')] 

in Jo 

= p(o) + i f dt'[v{t'),pm 



iti Jo 

1 



dt' / dt"[V(i'),[V(t"),p(0)]] + --- (8) 



(ifi)^ Ja ^0 
where 

p{t) EE e*^''*/^p(t)e-^^°*/^, V{t) EE e''^^ot/hy^-iHot/h^ 
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The reduced density matrix for the system oscillator is introduced as 

{ps)mn{t) = Tr{F„,„p(t)} = Tr{F,„„e-'«°*/''p(i)e*«°*/''}, (10) 

P„„ EE |n)(m|®lB, (11) 

p(0) = Ps(^PB^Y.^Ps)ki\k){l\(g>e-f'^^/ZB, (12) 

Zb = Tri3{e-'3^«}, (13) 

where the initial state is assumed to be a direct product state of ps and pB — e^'^^^/Zs, i.e., we have assumed 
that the bath is in thermal equilibrium. Here \k) is the vibrational eigenstate for the system Hamiltonian Hs, i.e., 
Ti.s\k) — Ek\k). If we assume that 7i/ is small, we can calculate \k) and Ek using the time-independent perturbation 
theory as shown in Appendix IXI 
We note that 

{ps),nn{t) = Tr{P,„„e-^^°*/''p(Oe^^°*/''} = TrB{p„„(i)}e-^"""*. (14) 
The lowest (second) order result for the density matrix is 

{ps)mn(.t) {ps)ilit) + (Ps)^'2(i) + iPsYrlit) + ' ' ' , (15) 

iPs)^°lit) = TrB{p,„„(0)}e— = (P5)™ne— (16) 
(P5)^LW - i f dt'TTB{(mmnpm\n)}e-^--' 
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= 77^ Z''^^' /' di"TrB{(m|[V(i'),[V(t"),p(0)]]|n)}e— 
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^' dt' ^' di" { (V™, {t')VM {t")) (ps);ne'("'"^*'+"-*") } e"^"-* 



1 



(zn)2 



..... e 
dt' /* dt" V f (Vi„(t')V™fe(t'')>(Ps)feie^'-'"''-*''+"'"*')| e-*"-* (18) 



where 



{Vkl{t)V„^n{t')) = TrB{psVfei(i)Km(t')}, (19) 
Vki{t) = (fc|V(OI0 = (fc|e^^«*/'"'Ve-*^«*/'^|0, (20) 
Uki = {Ek-Ei)/h. (21) 

Note that, in the above formulas, the time dependence is only induced by the bath Hamiltonian Tis- 
For the matrix elements of the interaction Hamiltonian V, we have 

{Vkl{t)) = ~{qs)kl{5T(t)) + {ql)u{5Q{t)), (22) 
{Vkl{t)Vn^n{t')) - {qs)kl{qs),nn{5T{t)5T{t')) + {ql)kl{ql),nn{5g{t)5g{t')) 

'[{qs)kl{ql),nn + {ql)kl{qs).nn]{5T{t)5g{t')) (23) 

where the value of {qs)ki and {q\)ki are given in Eqs. (jA16p - ljA2ip for the case of a cubic oscillator. Since {5J-) — 
and {5g) = 0, we have {Vki{t)) = and {psirkiit) = 0. 
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C. VER formula 



We first calculate the diagonal elements of the density matrix {ps)ii{t) {i — 0, 1) by assuming that the initial state 
is the first vihrationally excited state: ps = This is a typical situation for VER though VER from highly 

excited states can be considered The density matrix (p5)oo(0 is written as 



(p5)oo(i) - ^( dt' dt"Rc{(Vio(t')Voi(t")>e'""^*''*"^} 



where ujs is the anharmonicity-corrected system frequency given by Eq. (|A13|) . 
From Eq. ((23, we have 



(24) 



(ps)oo(i) - ^fe)?o / dt' 
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t' 



dt' I dt"Re\(Sg{t' -t")sg{o))e"^^'^*'-*"^\ 

Jo ^ ' 



(«s)io(g|)io f dt' r dt"Kc[{5T{t' -t")5g{0))e''''^''-'"^] 



2 

4 

" Jo Jo 

Using the explicit expressions for the correlation functions 7], the final VER formula is obtained as 

2 ^ 



(25) 



(26) 



where Utiyt) is defined as 



* , I , /I ^, I 1 — cosi7t 

dt' / dt"cosn{t' - t") = 

Jo 

t' 



Jo 



vt{n) = dt' dt" sin n{t' - t") = 



fit — sin fit 



(27) 
(28) 



and vt{^l) is defined for later use. The coefficients are defined in Appendix IbI Equation H26() is our final formula for 
VER. 

If we take the long time limit of this formula (which is equivalent to the Markov approximation), we obtain a 
formula for the VER rate 



^0^1 = ^(ps)oo(i) 



Q,/3 



27r 



(29) 



where we have used 



TTS{n). 



(30) 



If qs — qs and ujs = ujs, i-e., {J-') = (^J) = and / = 0, we recover the Maradudin-Fein formula from Eq. It 
follows that Eq. ^ is a generalization of the Maradudin-Fein formula, which can describe the time development of 
a density matrix. 



5 



D. Dephasing formula 

We calculate the off diagonal elements of the density matrix {ps)io{t) by assuming that the initial state is the 
superposition state between |0) and |1): ps = (1/2)(|0)(0| + |0)(1| + |1)(0| + |1)(1|) O. That is, (ps)fei = 1/2 for all 
k and /. This is a simplified situation to consider dephasing in a two- level system. We have 

ipsUit) - {ps)[l\t) + {ps)[l\t) + {ps)[l\t) + --- 

= le-''^-^*(l + r(i)(t)+r(2)(t) + ...). (31) 

By the definition of the interaction Hamiltonian f Appendix we have r'-^\t) = 0. The remaining term A^^{t) is 
decomposed as 

r(^)(t) = -rPAt)-r^GGit)~rPait), (32) 
rfpii) = ^[^^' [ di"Re{(Vio(i')Voi(t")>}[e^""^*'"*"^-e'""^*'+'"^] 



JO 



2 



- ft2fe)?oy dt' j di"Re{((5.F(<')^.F(i"))}[e^'^"(*'-*"^ -e^''^(*'+*")] 
+ 4('?s)?o f dt' f dt"Rc{((5g(t')<5^(t"))}[e'^"^*'-*")-e^^-'(*'+*")] 

- ^ (qs) loiql) 10 J^t' dt"Re{{dT{t')6g{t"))}[e''''^''-'"^ - e'^'^''+''\ (33) 
4gW = ^jy^' I ^i"{([Vii(0-Voo(t')]Vii(<")> + (Voo(t")[Voo(t')-Vii(t')]>} 

= ^[(95)11 -fe)oo]' fdt' f dt"Re{5Ht')5:F{t")) 



Jo 



^[(91)11 -(9|)oo]'y^ dt' dt"Re{6git')Sgit")) 
^[(95)11 -fe)oo][(g|)ii-(g|)oo] jy*' dt"Re{ST{t')Sg{t")) 
^[{ls)l, - fe)oo]^'rf^'/ dt"l^{5Ht')6Ht")) 
- (ls)lo] l^t' 1^ dt"Im{5g{t')5g{t")) 

2i 



^2[('Z5)ii(<z|)ii~(g|)oo(9|)oo] dt' dt"lm{ST{t')Sg{t")), (34) 

= 7l^[fe)ii-fe)oo]fe)io Trfi' r dt"lm{{dT{t')ST{t"))}[^'^-'' - e^''-'"] 
C^f^) Jo Jo 

t t' 

+ (^[(9s)ii-(95)oo](g^)io^ dt' dt"luv{{5g{t')5g{t"))W^^'' -e^^^'"] 

- jT^{[{qs)ii - (gs)oo](g|)io + - {ql)oo]{qs)w} 

X f dt' f dt"lra{{5T{t')6g(t"))}[e''^''' -e'^^'"] (35) 
Jo Jo 

where the subscript denotes that, e.g., for r^p(t), the dominant contribution comes from {5T{t)5T{Q)) when the 
effects of the bath and the system anharmonicity are both weak. 



6 



After similar calculations as done for the VER formula above, we obtain 



where 



(2) 



(2) 
^FG 



a, 13 



Q,/3 

i 



J2 [(D-f - i^+f )«t(c^a + u;p) + Di^J'vtiLO^ - up) + {Di^ - Dl^)vt{uj^) 



a, 13 



(^) = ^ E [(^-- - Ef+)9t{cOa + iop) + Ef_gt{LU^ - cop) + [El - E^)g,M 



a, 13 



dt' / dt"cosn{t' - t")[e'"^«(*'-*") - e»'^s(*'+*")] 
Jo 



1 



1 



Cjs — ^ (i)5 + r2j l_ (1)5 + ri 

5t(r2) = i / di' / dt" sinVl{t' - t")[e 



(Ds — f2 



2(1)5 



^0 



(I's(l + e*"^*)(l - cosm) - - e'"s*) sinm 
17((D| - f72) 



(36) 



(37) 
(38) 



(39) 



(40) 



and the coefficients are defined in Appendix IbI Equation (|31|l and Eqs. (|36|) - (|38|l constitute the dephasing formula. 

Dephasing properties are characterized by the decaying behavior of this off diagonal density matrix. Incidentally, 
as an alternative approach, one might use the von Neuman entropy or linear entropy for the reduced system as an 
indicator of dephasing [l^ . 



E. Frequency autocorrelation function 



Using the time-independent perturbation theory for the interaction, we obtain the fluctuation of the system fre- 
quency as 



SLu{t) = 



Vii(t)- Voo(i) _ {qs)ii-{qs) 



00 



dJ^it) 



(gi)ii - (g|) 



00 



h 



sg{t). 



Hence we have 



Re{Sijit')Su;in) = - ([(55)11 - (q s) oo? Re {mtWn) + K^Dii - (g|)oo]'Re(5t;(t')5t/(t")) 
-2[(q5)ii - (gs)oo][(9|)ii - {ql)oo]Re{STit')5git"))) . 

This turns out to be the second derivative of Re{r^^(^)}, i.e., 

C{t) EE Re{6uj{t)Scu{0)) = ^Re{r^^l{t)}. 
From this correlation function, we can define a pure dephasing time T2 as 

-=r 

Note that this is different from a correlation time defined by 

1 



C{t)dt. 



C{t)dt 



CiO) JO 

which leads to the relation l/Tj^* — C{0)tc. Note that Tj* ^^^d Tc are inversely related. 



(41) 



(42) 



(43) 



(44) 



(45) 



7 



III. APPLICATION: NMA IN HEAVY WATER 
A. NMA-D in heavy water 

We now apply our formulas to VER and dephasing problems of TV-methylacetamide (NMA) in heavy water. In 
many theoretical and experimental studies, this molecule CH3-NH-CO-CH3 is taken to be a model "minimal" peptide 
system because it contains a peptide bond (-NH-CO-). For example, Gerber and coworkers calculated the vibrational 
frequencies for this molecule using the vibrational self-consistent field (VSCF) method Nguyen and Stock worked 
to characterize VER in this molecule using a quasi-classical method Skinner and coworkers investigated the 

dephasing properties of the amide I mode using their correlation method combined with ab initio (DFT) calculations 
[T7|. Employing 2D-IR spectroscopy, Zanni et al. measured Ti and T2* for the amide I mode in this molecule, 
which were reported to be Ti ~ 0.45 ps and Tj* ~ 1.12 ps, whereas Woutersen et al. obtained T2* ~ 0.8 ps. 

In our study, we deuterate the system to NMAD/D2O so that the amide I mode, localized around the CO bond, 
can be clearly recognized as a single peak in the spectrum. In the following numerical calculations based on the 
CHARMM force field ^20|, its frequency is ~ 1690 cm~^ which fluctuates depending on the structure (see Fig. 
whereas the experimental and DFT values are 1717 cm~^ and 1738 cm~^, respectively [l7^ . 

B. Procedure 

We have applied the following general procedure. (1) Run an equilibrium simulation. (2) Sample several trajectories 
during the run. (3) Delete atoms of each configuration except the "active" region around the system oscillator. We 
introduce a cut off radius Rc around a certain atom within the active site. (4) Calculate instantaneous normal modes 
(INMs) [231 for each "reduced" configuration, ignoring all the imaginary frequencies (5) Calculate anharmonic 
coupling elements with the finite difference method 7] using the obtained INMs. (6) Insert the results in the VER 
formula [Eq. if^ ] and dephasing formula [Eq. ll^ with Eqs. (7) Ensemble average the resultant density 

matrix, and estimate the VER and dephasing rates (times), if possible. 

This procedure seems to be straighforward. However, when applied to real systems like peptides or proteins, we need 
to carefully treat the the effects of the bath. In Fig. ^ we plot the system frequency uos for 100 sample trajectories 
from an equilibrium run at 300K. We see that the amide I mode frequency changes depending on the structure; the 
deviation can amount to 1%. 

1698 I 1 1 1 1 1 1 1 1 1 1 




10 20 30 40 50 60 70 80 90 100 

sample index 



FIG. 1: Instantaneous normal mode frequency of the system uos for 100 different sample trajectories at 300K where the cut 
off radius is Rc — 10 A. 

Furthermore the frequency is renormalized according to Eqs. © and (|A13|I . and such an effect can be anomalously 
large if we include all the contribution from low frequency components. Hence we need to introduce a cutoff frequency 
Wc, below which the contribution is neglected. This is physically sound, because we are dealing with time-dependent 
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phenomena, and such low frequency components correspond to longer time behavior. However, we are now interested 
in rather short time dynamics, so such contributions should not play a role. In fact, the final result of VER does 
not depend much on the choice of uJc, whereas that of dephasing does. We need to admit that for now this is just a 
remedy. We discuss how to improve this situation later. 

C. VER properties of NMA-D 

First we consider the VER properties of the amide I mode as shown in Fig. |21 We use the following relation 

piiit) = 1 - poo(t) = cxp[-s(i)] (46) 

and hypothesize that s£i) ~ Poo(0: which is definitely true when poo{t) <^ 1, and might be justified using the cumulant 
expansion technique [23, ll^ . 

We calculated the density matrix for the following three cases: (a) NMA-D in heavy water with CHARMM force 
field at 300K, (b) NMA-D in vacuum with CHARMM force field at OK, and (c) NMA-D in vacuum with DFT 
force field at OK. Here we have used Rc — 10 A and ujc = 10 cm^^ for case (a). The result for VER does not 
depend sensitively on these parameters. For cases (b) and (c), we must take a special care: It is known that the low 
frequency components cause serious problems for vibrational frequency calculations fSJ , so we need to eliminate the 
low frequency components. In this work, we exclude several normal modes if their frequency is less than 300 cm~^. 
See Table n 

TABLE I: Normal mode frequencies (in cm"^) for ah initio (left) and CHARMM (right) NMA. The level of the ab initio 
calculation is B3LYP/6-3H-G(d). 



Mode index a 


UJc (ab initio) 


LOc (CHARMM) 


Mode index a 


uja (ab initio) 


UJc, (CHARMM) 


1 


31.5 


64.1 


16 


1417.9 


1380.2 


2 


71.6 


88.9 


17 


1436.1 


1408.7 


3 


170.2 


192.3 


18 


1483.6 


1415.4 


4 


259.6 


269.5 


19 


1495.1 


1418.5 


5 


282.3 


426.7 


20 


1499.4 


1425.7 


6 


421.9 


536.3 


21 


1516.7 


1444.7 


7 


619.1 


575.9 


22 


1535.9 


1563.1 


8 


619.9 


741.9 


23 


1745.9 


1678.1 


9 


868.7 


757.0 


24 


2671.1 


2445.0 


10 


946.1 


854.4 


25 


3058.5 


2852.8 


11 


1012.4 


964.3 


26 


3058.9 


2914.3 


12 


1066.1 


1055.9 


27 


3116.5 


2914.8 


13 


1144.5 


1075.7 


28 


3130.9 


2917.2 


14 


1158.9 


1088.5 


29 


3135.9 


2975.3 


15 


1207.6 


1123.5 


30 


3148.9 


2975.5 



By using a fitting form s{t) — t/Ti, where Ti is the VER time, we estimate that Ti ~ 0.5 ps at 300K and 0.6 
ps at OK from the initial decay. The former estimate is rather similar to the experimental value Ti ~ 0.45 ps 
whereas Nguyen-Stock's quasi-classical estimate is Ti ~ 1.5 ps fT6j. Considering that the estimate at 300K is rather 
close to that at OK, we can conclude that quantum effects are important to describe VER for the amide I mode of 
NMA-D in heavy water. However, the decay at the later stage becomes very slow at OK as expected because there is 
no environment. (In the vacuum cases, we only use one minimized structure, thus there is no ensemble average, and 
the oscillatory behavior remains.) 

The results are similar for NMA-D in vacuum with different force fields. It is known that NMA with CHARMM 
force field is not well characterized around the methyl groups , but this fact does not affect the VER properties of 
the amide I mode. 

We have analyzed the mechanism of VER in terms of the VER pathway. In Table m we show several mode 
combinations that contribute most to s{t) for NMA-D in heavy water at 300K. These eigenvectors (normal modes) 
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are well localized around NMA (Fig. O, especially on the CO bond (Table ^HJ. There is very little contribution from 
the surrounding water. (This is expected from the previous result of Kidera and coworkers [26] ■) Similar "resonant" 
mode combinations can be found in the isolated NMA-D cases. See Tables Hvl and Ivl This means that the initial 
stage of VER of NMA-D in heavy water is dominated by intramolecular vibrational redistribution (IVR) localized 
near the peptide bond. This result might explain why the amide I mode, in many peptide systems with differing 
environments, appears to have similar VER times |27l |. Note that this is the case for a localized mode such as the 
deuterated amide I mode. A collective mode can decay with a different VER pathway, as shown by Austin's group 

M. 




0.7 - 
0.6 - 
0.5 
0.4 
0.3 
0.2 



T 1 1 1 1 r 

DFT(OK) — ^ 
CHARMM(OK) — x- 
CHARMM (300K) i b 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

?(ps) 



FIG. 2: Time evolution of the excited density matrix for NMA-D in vacuum (DFT and CHARMM) at OK, and for NMA-D 
in heavy water (CHARMM) at 300K. The level of DFT is B3LYP/6-31+G(d). 



TABLE IL The most dominant VER pathways for the amide I mode of NMA-D in heavy water. Act) is defined by \ujs—^a —^p 



Mode combination (q, fS) 


frequency (cm ^) 


Contribution to s(f) 


Atj (cm~^) 


1143 + 1143 


778.6 4- 778.6 


0.04 


125.7 


1147 + 1134 


1085.3 570.0 


0.04 


27.6 


1147 + 1135 


1085.3 -f 570.9 


0.01 


26.6 


1147 + 1136 


1085.3 578.8 


0.02 


18.8 


1147 + 1137 


1085.3 -^- 581.3 


0.03 


16.2 


1147 + 1140 


1085.3 612.1 


0.46 


14.5 


1148 + 1132 


1127.8 + 558.5 


0.01 


3.4 


1148 + 1134 


1127.8 -1- 570.0 


0.11 


14.9 


1148 + 1135 


1127.8 570.9 


0.03 


15.9 



D. Dephasing properties of NMA-D 

We now consider the dephasing properties of the amide 1 mode. The off-diagonal density matrix is written as 

Pio(t) - - rfpify - r^^lit) - r%{t)] ~ ^ e-'^-*-'---'*)-''^'^'*)-''-'^^*) (47) 

and we analyze each contribution to the density matrix seperately. 
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T ffr 



^ 0.2 - 
0.1 - 



I 1 1 ^ ++ 1^ 1 ^v:' 1 1 J 

-500 500 1000 1500 2000 2500 3000 3500 

(Oa {cm ') 



FIG. 3: Norm of the eigenvectors (normal modes) with the exception of the contribution from NMA-D, which is defined by 
'}2i£Wa.tcii-'^i + + where i comes from water degrees of freedom alone. 



TABLE III: The most localized modes around the CO bond in NMA-D. The norm is defined by I]igcobond(^i + Vi + ^?)- 
The 1345th mode is the amide I mode. 



Mode index a 


frequency (cm ^) 


Contribution to norm 


1134 


570.0 


0.14 


1140 


612.1 


0.26 


1142 


771.2 


0.35 


1143 


778.6 


0.41 


1146 


1013.6 


0.13 


1147 


1085.3 


0.26 


1148 


1127.8 


0.17 


1340 


1452.1 


0.36 


1345 


1689.6 


0.92 



TABLE IV: The most dominant VER pathways for NMA-D in vacuum with the ab initio potential (B3LYP/6-31+G(d)). 



Mode combination (q, I3) 


frequency (cm ^) 


Contribution to s{t) 


Alo (cm~^) 


9 + 9 


868.7 + 868.7 


0.13 


2.7 


13 + 8 


1144.5+ 620.0 


0.02 


29.8 



TABLE V: The most dominant VER pathways for NMA-D in vacuum with the CHARMM force field. 



Mode combination (a, l3) 


frequency (cm ^) 


Contribution to s{t) 


Auj (cm ^) 


8 + 8 


741.9 +741.9 


0.02 


194.1 


14 + 6 


1088.5+536.3 


0.08 


53.1 


14 + 8 


1088.5+741.9 


0.02 


152.5 


15 + 7 


1123.5+575.9 


0.08 


21.6 
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In Fig. 21 we show the result with i?c = 10 A and uoc = 10 cm ^. We can see that the following relation holds 

Re{^FfW}-sW/2-V(2ri). (48) 
If we further assume that Re{r[?^(t)} ~ t/TJ^* and Re{r^Q(i)} ~ 0, we have 

^^^ + ^ (49) 
T2 2Ti T*' ^ ' 

This is a standard expression connecting Ti and T2 |29j | , and holds under the Markov approximation. We can see that 

(2) (2) 

I{je{rpQ{t)} ~ holds, but it is difficult to judge whether Re{rQQ(t)} ~ i/Tj* holds or not. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
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FIG. 4: Dephasing properties of the amide I mode of NMA-D in heavy water. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

tips) 



0.045 

0.04 
0.035 

0.03 
0.025 

0.02 
0.015 

0.01 
0.005 


-0.005 





(i), = 1.0 cm"' 






0)^=10.0 cm-' 


---X--- - 




ffl,= 100.0 cm-' 
















\\ 






~\ \ 

























0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

tips) 



FIG. 5: Left; VER (left) and frequency auto correlation calculations (right) at 300K with different cut-off frequencies ujc- 

There are more serious problems: as mentioned by Mikami and Okazaki |lO|. the diagonal terms contribute most 
for dephasing, i.e., the second term in Eq. H37() is a dominant contribution for dephasing. Furthermore, the coefficients 
are dominant factors, so this means that the low frequency (and thus delocalized) modes contribute most. In this 
paper, we have employed two cut-off parameters: Rc and lUc- If Rc is large enough, it is fine, but the choice of uJc can 
be arbitrary. Figure [S] shows the dependence of the results on uic- The VER results do not depend on the choice of 
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Wc because there is a resonant condition which should be met, but the dephasing results do. We need to be cautious 
in the interpretation of our results for dephasing. One way to get rid of this problem is to go back to the original 
expression Eq. using the force and force-constant autocorrelation functions {6!F{t)6!F{0)) and {SG{t)SG{0)) . Here 

(2) 

^gg(^) calculated as time integral of these correlation functions, which can be calculated using classical mechanics. 
This is in the same spirit as the quantum correction factor method |30| |. which is an approximation to quantum effects. 
In this case, we only need to consider the zero frequency component, so the classical mechanics should work well and 
quantum effects should be less important. 

E. Discussions 

We found that there are several resonant modes in NMA-D, which form the main VER pathways within the 
molecule. Gerber and coworkers reported that the amide I mode in NMA is very weakly coupled to other modes . 
We expect that this discrepancy results from (a) the use of only pair interactions between normal modes to reduce the 
computational cost, (b) the level of the ab initio method: they used MP2/DZP whereas we used B3LYP/6-31+G(d), 
and (c) the criterion of the mode-mode coupling: their criterion is not directly related to VER. 

It is important and interesting to clarify the nature of VER in the amide I mode in more detail. Note the importance 
of the system anharmonicity. The effect of the system anharmonicity defined in Eq. ljA10|) is very weak for the 
CHARMM case: e = 10~^, but it is not for the ab initio case: e = 10~^. According to Eq. I) A 13(1 . this anharmonicity 
shifts the system frequency by 0.6%, which amounts to 10 cm"-'^. The resonant condition changes compared to the 
case without anharmonicity. Of course, dephasing is also affected by this amount of anharmonicity. To address these 
issues, we must develop QM/MM type methods, which will be described elsewhere. Another interesting system to 
investigate anharmonicity is a highly excited bond such as the highly excited CO bond in myoglobin [T^ . 

It is important to assess our strategy: our perturbative expansion and cut-off strategy are approximate. It would 
be profitable and interes ting to compare this strategy with others, including the time-dependent vibrational self- 
consistent field methods [3in3^. the semiclassical method js^], and the path integral method [l^- The application 
of our strategy to protein systems, including cytochrome c, will be described elsewhere [35|. 

IV. SUMMARY 

In this paper, we have derived formulas for VER and dephasing for an anharmonic (cubic) oscillator coupled to a 
harmonic bath through 3rd and 4th order coupling elements. We employed time-dependent perturbation theory and 
did not take the infinite time limit as is done in the derivation of the Maradudin-Fein formula. Hence our formulas 
do not assume the Markov properties of the system, and can describe short time behavior that can be important for 
VER and dephasing properties of localized modes in peptides or proteins. Our final results are the VER formula 
[Eq. H2()|l ] and dephasing formula [Eq. H31I) with Eqs. I|36|) - (|38ll ]. As a test case, we have studied the amide I mode of 
iV-methylacetamide in heavy water. We found that the VER time is 0.5 ps at 300K, which is in good accord with the 
experimental value, and clarified that the VER mechanism is mainly localized around the peptide bond in NMA-D; 
VER is dominated by IVR within the molecule. We also investigated the dephasing properties of the amide I mode, 
and met with some problems. We proposed a new method to overcome these problems using classical correlation 
function calculations. 
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APPENDIX A: SYSTEM PARAMETERS FOR A CUBIC OSCILLATOR 

We assume that the system-bath interaction can be Taylor expanded using the bath coordinate qo„ and that the 
fluctuating force and the fluctuating force constant can be expressed as 

= ^Csa/3(gag/3 - (feg/j)), (Al) 
SG = ^Cssa/siqaQp - {qaQp)) + ^Cssaqa- (A2) 

a,f3 a 

In real molecular systems such as peptides or proteins, the coefficients in V and the anharmonicity parameter in 
Hf are calculated as 

1 d^V 

Csap = --^-^ — ^ — — , (A3) 

2 dqsoqadqp 

1 d'^V 



1 d'^V 

CsSal3 = T ^ 2 fl a ' (^^) 

4 dq^dqadqn 

dqi 



where V represents a potential function for the system considered. This potential function can be an empirical force 
field (CHARMM, Amber) or an ab initio potential calculated by any level of theory. 

Assuming that the cubic anhamonicity / in the system is small, we use the time-independent perturbation theory 
to calculate the eigen energies and vectors. We quote from J.J. Sakurai ^3^: 




p(0) p(0) 



T4.V^„ ^i^(o). I (A8) 



where En^ = hll)s{n + 1/2), jfc^"^) is the k-th eigenfunction of the harmonic oscillator, and 



hioss 



n{n - l)(n - 2)4, n-3 + 3ny/n6, 



k,n-- 1 



+3(n + l)Vn+lSk,„+i + ^/in+^){n + 2){n + 3)Sk,n+3\ (A9) 

where 

1 f f ^ \ ^^'^ 
tiws 6 \2ujs ) ^ 

is a dimensionless paramater representing the strength of the anharmonicity of the system. Note that Vkn becomes 
nonzero only when |fc — n| = 1 or |/c — n| = 3. 
We explicitly have 

^^^fi^_r^_^^^^^_ 22^2), (All) 

2 hujs Shujs 2 

^ 3^ ^ ^ _ ^ _ 1^ ^ 142.^,. (A12) 

2 Thus n^s Anios 2 



14 



The anharmonicity-corrected frequency is 

us=^^^=u>sil-60s^). (A13) 
Next we calculate the matrix elements for qs and qg. We write the eigenfunctions: 

fe=i,3 -C'o -^k k,ieSo K^o - Jy'^o ~ ) 

ID = 11(0.) + y: i'''°'>:^#>t+ e i'''°'> ,^„i ^.t";.,, ^,0,, (A15) 

fe=o,2,4 ^1 --f^fe k,ieSi - ^k - ) 

where So represents (/ = 1, A: = 2) or = 1, A: = 4) or (/ = 3, /c = 2) or (/ = 3, /c = 4) or (/ = 3, A: = 6), and 5*1 does 
(; = 0, fc = 3) or = 2, fc = 3) or (/ = 2, A: = 5) or = 4, A: = 3) or (I = 4, A; = 5), or (I = 4, A; = 7). Note that these 
eigenvectors are not normalized, so we need to renormalize them before or after calculations. 
After some lengthy but straighforward calculations, we have 

fe)io = (l|gs + &|0) = (gs)oi =a(l + 22e2), (A16) 

fe)oo = ms + b\0) = h-Qae, (A17) 

(95)11 = {l\qs + h\l) = b-l^ae, (A18) 

(«|)io = (I|fe + 6)'|0) = ((7|)oi=2a6-20a2£ + 44a6e^ (A19) 

(al)oo = (0|(gs + b)^|0) = a2 + 62-12a6£ + 88aV, (A20) 

= (l|(gs + 6)2|l) =3a2 + 62-36a6e+568aV (A21) 

where 



a = ,/^ (A22) 



is the fundamental length charactering the system oscillator. 
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APPENDIX B: THE COEFFICIENTS USED IN THE FORMULAS 



Using the expression derived previously for the force-force correlation function , the coeiEcients in our VER and 
dephasing formulas are expressed as 

C"'^ = (^^3 ^S) ={fe)ioC^^"/3-('?l)ioCsso/5}'s"^ (Bl) 
= {[{qshi - {qsMCsap - [(4)ii - (g^)oo]Cssa;3}'s"^ (B2) 

= {[(9s)ll - {qs)0Q]CsaP - [(<7s)ll - {qs)Q0]CsSap} 

X {[fe)ii + iqs)oo]Cso.0 - [(g|)ii + {ql)oo]Cssa0} S"^ (B3) 



E 



a/3 



E^^_ Ef_ 



{[(9s)ll - {qs)0Q]CsaP - [{ql)ll - {qs)oo]CsSap} 

X {{qsUCscp - {qDioCsSc^} S"^ (B4) 

ft' f {I + na,)il + np) 2{1 + no,)nf3 \ ^ ^^^^ 

D""" - I 1 = - (9s)oo]'C|5„R", (B7) 



jQ/3 _ 



C" = I 1 = (9|)?oC|saR", (B6) 



D'" = t^'^^^" " (95)00] + (g|)oo]C|5„R", (B8) 

E" = ( ) = " (9l)oo](9^)ioC|s„R", (B9) 



2ti)Q \ rio 

where Ua = l/(e''''"° — 1) is the thermal phonon number. 

To calculate Cjs and h in Eqs. © and lO, we use the following 



{Ht)) = im) = f E + 2"-)' (Bii) 

a 
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